• TIDYVERSE(DPLYR):
  • Manipulación de datos. Similitud a PANDAS (Python)
  • Flujo de trabajo
  • Permite mayor simplicidad en la codificación
  • Concatenación (pipe)
  • Sintaxis entendible
  • Similar a manipulación vía SQL
  • Desventaja: mucho uso de memoria.
• GGPLOT:
  • Gramática de los gráficos por capas (http://vita.had.co.nz/papers/layered-grammar.pdf)
  • Permite seguir los pasos y los distintos componentes del gráfico
  • Permite una fácil forma de iteración para distintos gráficos

Procesos: bajar base de datos (1), seleccionar columnas (dplyr::select) (2), formatear fecha (desde formato de texto día-mes-año a “Ymd”) (3), filtrar regiones (sólo Coquimbo y Valparaíso) (4), generar figura de serie de datos (Fecha como eje X y Nuevo Caso Confirmado como eje Y), discriminando por Región por colores (5), generamos un área gris del rango intercuartil del resto de regiones, excluyendo a la región metropolitana (al 80% o fill=grey70) con transparencia para que no se superponga (alpha=.4) (6), definimos los colores de las líneas: rojo para Coquimbo, azul para Valparaíso y negro para la mediana del resto de regiones (scale_color_manual) (7).

covid19_chile_coq_val <- 
    readr::read_csv("https://raw.githubusercontent.com/ivanMSC/COVID19_Chile/master/covid19_chile.csv") %>% 
    dplyr::select(Fecha, Region, `Nuevo Confirmado`) %>%
    dplyr::mutate(Fecha=as.Date.character(Fecha,tryFormats = "%d-%m-%Y")) %>% 
    dplyr::filter(Region %in% c("Valparaíso","Coquimbo")) %>% 
    ggplot()+
      geom_line(aes(x=as.Date(Fecha), y=`Nuevo Confirmado`, color=Region))+
      scale_x_date()+
      xlab("Fecha")+
      sjPlot::theme_blank()+
      theme(legend.position="bottom")

## 
## -- Column specification --------------------------------------------------------
## cols(
##   Fecha = col_character(),
##   Region = col_character(),
##   `Nuevo Confirmado` = col_double(),
##   `Nuevo Muerte` = col_double(),
##   `Nuevo Recuperado` = col_double(),
##   `Acum Confirmado` = col_double(),
##   `Acum Muerte` = col_double(),
##   `Acum Recuperado` = col_double(),
##   `Casos UCI` = col_double()
## )

covid19_chile_resto <-
    readr::read_csv("https://raw.githubusercontent.com/ivanMSC/COVID19_Chile/master/covid19_chile.csv") %>% 
    dplyr::select(Fecha, Region, `Nuevo Confirmado`) %>%
    dplyr::mutate(Fecha=as.Date.character(Fecha,tryFormats = "%d-%m-%Y")) %>% 
    dplyr::filter(!Region %in% c("Valparaíso","Coquimbo")) %>% 
    dplyr::filter(!Region %in% c("Metropolitana")) %>% 
    dplyr::group_by(Fecha)%>% 
    dplyr::summarise(mdn=median(`Nuevo Confirmado`,na.rm=T),
                     q1=quantile(`Nuevo Confirmado`,.25,na.rm=T),
                     q3=quantile(`Nuevo Confirmado`,.75,na.rm=T))

## 
## -- Column specification --------------------------------------------------------
## cols(
##   Fecha = col_character(),
##   Region = col_character(),
##   `Nuevo Confirmado` = col_double(),
##   `Nuevo Muerte` = col_double(),
##   `Nuevo Recuperado` = col_double(),
##   `Acum Confirmado` = col_double(),
##   `Acum Muerte` = col_double(),
##   `Acum Recuperado` = col_double(),
##   `Casos UCI` = col_double()
## )

covid19_chile_coq_val+ 
    geom_ribbon(data= covid19_chile_resto, aes(y= mdn, x= Fecha, ymin = q1, ymax = q3), fill = "grey70", color="grey70", alpha=.4)+
    geom_line(data= covid19_chile_resto, aes(x= Fecha, y= mdn, color="black"))+
    scale_color_manual(name= "Región", values = c("black","red","blue"), 
                       label=c("Mediana (resto regiones\nexcluyendo Metropolitana)",
                               "Coquimbo",
                               "Valparaíso"))

• ¿No se ve más ordenado?, lo valorarán cuando trabajen con condiciones anidadas (IFs concatenados)

Factorizamos lo que tienen en común las bases de datos, y para cada capa filtramos los datos de interés.

covid19_chile <-
    readr::read_csv("https://raw.githubusercontent.com/ivanMSC/COVID19_Chile/master/covid19_chile.csv") %>% 
    dplyr::select(Fecha, Region, `Nuevo Confirmado`) %>%
    dplyr::mutate(Fecha=as.Date.character(Fecha,tryFormats = "%d-%m-%Y"))

## 
## -- Column specification --------------------------------------------------------
## cols(
##   Fecha = col_character(),
##   Region = col_character(),
##   `Nuevo Confirmado` = col_double(),
##   `Nuevo Muerte` = col_double(),
##   `Nuevo Recuperado` = col_double(),
##   `Acum Confirmado` = col_double(),
##   `Acum Muerte` = col_double(),
##   `Acum Recuperado` = col_double(),
##   `Casos UCI` = col_double()
## )

ggplot()+
      geom_line(data=covid19_chile %>% 
                  dplyr::filter(Region %in% c("Coquimbo","Valparaíso")),
                aes(x=as.Date(Fecha), y=`Nuevo Confirmado`, color=Region))+
      scale_x_date()+
      xlab("Fecha")+
      sjPlot::theme_blank()+
      theme(legend.position="bottom")+
      geom_ribbon(data= covid19_chile %>% 
                  dplyr::filter(!Region %in% c("Coquimbo","Valparaíso","Metropolitana")) %>% 
                  dplyr::group_by(Fecha)%>% 
                  dplyr::summarise(mdn=median(`Nuevo Confirmado`,na.rm=T),
                     q1=quantile(`Nuevo Confirmado`,.25,na.rm=T),
                     q3=quantile(`Nuevo Confirmado`,.75,na.rm=T)), 
                aes(y= mdn, x= Fecha, ymin = q1, ymax = q3), fill = "grey70", color="grey70", alpha=.4)+
    geom_line(data= covid19_chile_resto, aes(x= Fecha, y= mdn,color="black"))+
    scale_color_manual(name= "Región", values = c("black","red","blue"), 
                       label=c("Mediana (resto regiones\nexcluyendo Metropolitana)",
                               "Coquimbo",
                               "Valparaíso"))

Muy complejo y todo, pero no entendí nada los gráficos anteriores. Si queremos separar los datos por región:

ggplot()+
  # Defino la línea a aprtir de datos con información de los casos, dejando sólo 2 regiones
      geom_line(data=covid19_chile %>% 
                  dplyr::filter(Region %in% c("Coquimbo","Valparaíso")),
  # Defino cómo serán interpretados los datos, eje x es la fecha e y es el número de confirmados
                aes(x=as.Date(Fecha), y=`Nuevo Confirmado`, color=Region))+
  # Como comparación, dejé una cita/área gris con información, excluyendo a Coquimbo y Valparaíso, además de una región muy influyente: Metropolitana (por su población). Dejamos como fuente de datos por la mediana de esas Regiones por cada día
        geom_ribbon(data= covid19_chile %>% 
                  dplyr::filter(!Region %in% c("Coquimbo","Valparaíso","Metropolitana")) %>% 
                  dplyr::group_by(Fecha)%>%
                  dplyr::summarise(mdn=median(`Nuevo Confirmado`,na.rm=T),
                     q1=quantile(`Nuevo Confirmado`,.25,na.rm=T),
                     q3=quantile(`Nuevo Confirmado`,.75,na.rm=T)), 
    # Definimos los elementos estéticos del gráfico, en el que dejamos el límite inferior en el percentil 25 y el límite superior en el percentil 75, coloreando el área 70% gris, aunque con una transparencia del 40%.
                aes(y= mdn, x= Fecha, ymin = q1, ymax = q3), fill = "grey70", color="grey70", alpha=.4)+
  # El eje x estaba en formato de fecha
      scale_x_date()+
  # Se añade una etiqueta "Fecha"
      xlab("Fecha")+
  # Se define un tema estético para todo el gráfico
      sjPlot::theme_blank()+
  # Definimos dónde dejar la etiqueta
      theme(legend.position="bottom")+
  # Dividir el gráfico por región de interés (en este caso, toma como referencia la primera capa [línea])
      facet_wrap(~Region, ncol = 1)

• Existen antecedentes de que gran parte de las parejas en Coquimbo y Valparaíso asistió a pubs y restaurantes dado que estaban en mayor proporción abiertos. El resto, recurrió a otros medios de contacto
• Esto no ocurrió en regiones, quienes se encontraban en cuarentena estricta
• Dado que algunos expertos señalan que la transmisibilidad se reduce con la primera dosis, mayores vacunaciones debiesen confundir las asociaciones
• ¿El 14 de febrero en Coquimbo y Valparaíso tuvo un efecto en el aumento de casos (controlando por vacunaciones)?
• Debemos incorporar las vacunaciones a la base de datos de nuevos casos
• Otro paquete de la familia Tidyverse es “tidyr”
• Similar a las funciones utilizadas por reshape (cambiar el formato largo y ancho de los datos)
• También mostraremos como ejemplo los join.
• Permiten combinar bases de datos
  • Left
  • Inner
  • Right
  • Otros (ej., probabilístico)

• “scientific generalization is a broader question than mathematical description.”
• No importa tanto el estimador puntual, si no su varianza (intervalos de confianza, márgenes de error MEDIBLES, poder estadístico).
• Hay muchas fuentes de sesgo o error sistemático que hay que atender, e incluso en distintas etapas del proceso de investigación.
• Aproximación cuasi-experimental: Distinguir efecto causal y resultado bajo observación pasiva. No puedo aislar ni manipular variables (entre ellas, la aleatorización de unidades)
• En algunos casos se utiliza para predecir, pero es distinto poner énfasis en el predictor (B o beta, inferencia causal) que en el valor predicho (Y) (Ver diferencias entre asociación predivtiva y una relación causal).

Figura. Varianza vs. Sesgo (Fuente: Singh, 2018, Mayo 21)

Panel

• Datos de panel: unidad de observación repetida en el tiempo
• Balanceados vs. Desbalanceados
• Las unidades tienen patrones de comportamiento individuales

library(DiagrammeR) 
gr1<-
DiagrammeR::grViz("
digraph grafiquito  {

# Nodes
  node [shape = plaintext]
  a [label = 'Confirmados \nPre 14 Feb\nZ(0)',fontsize=10]
  b [label = 'Coquimbo y\nValparaíso\n(A)',fontsize=10]
  c [label = 'Confirmados \nPost 14 Feb\nZ(1)',fontsize=10]
  d [label = 'Confusores no-medidos\n(U)',fontsize=10]
  e [label = 'Confusores medidos\nPre 14 Feb\n(C)',fontsize=10]
  blank [label = '', width = 0.01, height = 0.01]
  g [label = '', fontsize=10, width = 0.001, height = 0.001, color=White]

# Edges
  edge [color = black]

  b -> c 
  b -> blank [rankdir = TD; arrowhead = none; color = white]
  blank -> d [arrowhead = none; color = white]
  a -> c [dir= both; style = dashed]
  d -> { b a c }
  e -> { b a c }
#  d -> g [ dir = none,  color = 'white',fontcolor = white,shape=none, width=0, height=0];

  { rank = same; rankdir = LR; b; a; c }
  { rankdir = TD; a; d; e }

# Graph
  graph [overlap = true]
}")
gr1

#Fig. 2. Directed acyclic graph depicting the causal association between the treatment A, pre-exposure outcome Z(0), post-exposure outcome Z(1), measured pre-exposure confounders C, and unmeasured confounders U.

• No podremos ver (U). Para ello debemos hacer pruebas de sensibilidad a confusores no-medidos

• Se estandarizaron las variables de acuerdo a la población de la región, por cada 100,000 hab.
• Se generó un modelo de efectos fijos para capturar eventuales cambios en coquimbo en comparación a otras regiones a partir del 14 de febrero
• ¿Fue la forma correcta de incorporar las bases de datos? ¿Se perdió información?

library(tidyr)
library(lmtest)
library(plm)
library(broom)
library(kableExtra)

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#Generación y compilación de las bases de datos
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

covid19_chile_vacuna_vs_confirmados <-
    readr::read_csv("https://raw.githubusercontent.com/MinCiencia/Datos-COVID19/master/output/producto76/vacunacion.csv") %>%
  #Formato largo de fechas
    tidyr::pivot_longer(cols=starts_with("20"), names_to="fecha", values_to="recuento") %>% 
  #Formato ancho por tipo de dosis (diferenciando primera y segunda en columnas distintas)
    tidyr::pivot_wider(names_from=starts_with("Dosis"), values_from="recuento") %>% 
    dplyr::mutate(fecha=as.Date.character(fecha,tryFormats = "%Y-%m-%d")) %>% 
  # Unir con la base de datos de casos por Región
    dplyr::left_join(covid19_chile, by=c("fecha"="Fecha", "Region"="Region")) %>% 
    dplyr::filter(Region!="Total") %>% 
    dplyr::select(-Segunda) %>% 
    dplyr::filter(!is.na(Primera)) %>% 
    dplyr::mutate(met=factor(ifelse(Region %in% c("Coquimbo", "Valparaíso"),1,0))) %>% 
    dplyr::mutate(txtime=factor(ifelse(fecha>="2021-02-14",1,0))) %>% 
    dplyr::rename("conf"="Nuevo Confirmado")

## 
## -- Column specification --------------------------------------------------------
## cols(
##   .default = col_double(),
##   Region = col_character(),
##   Dosis = col_character()
## )
## i Use `spec()` for the full column specifications.

#Base de datos de PCR's y UCI's
covid19_chile_vacuna_vs_confirmados %>%     
  dplyr::left_join(readr::read_csv("https://raw.githubusercontent.com/MinCiencia/Datos-COVID19/master/output/producto7/PCR_std.csv"), 
                   #Los criterios para parear es la Región(i) y la fecha(t) 
                   by=c("Region", "fecha")) %>% 
  dplyr::left_join(readr::read_csv("https://raw.githubusercontent.com/MinCiencia/Datos-COVID19/master/output/producto8/UCI_std.csv") %>% 
                     #sacar variables redundantes: código de región, población y número
                     dplyr::select(-`Codigo region`,-Poblacion), 
                   by=c("Region", "fecha")) %>% 
  dplyr::mutate(conf_100mil=(conf/Poblacion)*100000) %>% 
  dplyr::mutate(Primera_100mil= (Primera/Poblacion)*100000) %>% 
  dplyr::mutate(PCR_100mil= (numero.x/Poblacion)*100000) %>% 
  dplyr::mutate(Hosp_UCI_100mil= (numero.y/Poblacion)*100000) %>% 
  dplyr::mutate(did=factor(as.numeric(met)*as.numeric(txtime)),
                did=factor(ifelse(did==4,1,0))) %>% 
  as.data.frame() %>% 
  #Para definir la base de datos
  assign("covid19_chile_vacuna_vs_confirmados_std_pcr_uci",., envir=.GlobalEnv)

## 
## -- Column specification --------------------------------------------------------
## cols(
##   Region = col_character(),
##   `Codigo region` = col_character(),
##   Poblacion = col_double(),
##   fecha = col_date(format = ""),
##   numero = col_double()
## )

## 
## -- Column specification --------------------------------------------------------
## cols(
##   Region = col_character(),
##   `Codigo region` = col_double(),
##   Poblacion = col_double(),
##   fecha = col_date(format = ""),
##   numero = col_double()
## )

paste0("La fecha mínima pareada fue ",min(covid19_chile_vacuna_vs_confirmados_std_pcr_uci$fecha))

## [1] "La fecha mínima pareada fue 2020-12-24"

paste0("La máxima fecha disponible fue ",max(covid19_chile_vacuna_vs_confirmados_std_pcr_uci$fecha))    

## [1] "La máxima fecha disponible fue 2021-03-06"

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
# Modelo de efectos fijos
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

did.reg_basico <- plm(conf_100mil ~ txtime + met + txtime:met, 
               data = covid19_chile_vacuna_vs_confirmados_std_pcr_uci,  index=c("Region", "fecha"), model = "within")
did.reg <- plm(conf_100mil ~ txtime + met + txtime:met + Primera_100mil + PCR_100mil + Hosp_UCI_100mil, 
               data = covid19_chile_vacuna_vs_confirmados_std_pcr_uci,  index=c("Region", "fecha"), model = "within")
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

modelo_simple<-
parameters::model_parameters(did.reg_basico) %>% 
  data.frame() %>% 
  #Cambiamos los nombres de las variables
  dplyr::mutate(Parameter= dplyr::case_when(Parameter=="txtime1"~"14 de Febrero hasta el día de obtención",
                                       Parameter=="Primera_100mil"~"Primera dosis por cada 100 mil hab.",
                                       Parameter=="PCR_100mil"~"PCR's por cada 100 mil hab.",
                                       Parameter=="Hosp_UCI_100mil"~"Hospitalizaciones por cada 100 mil hab.",
                                       Parameter=="txtime1:met1"~"Coquimbo y Valparaíso x 14 de Febrero"
                                       )) %>% 
  dplyr::mutate_at(vars("Coefficient","SE","CI_low", "CI_high","t", "df_error"), funs(round(.,2))) %>% 
  #Aproximamos a 3, de manera que los interesados puedan ver el valor p con 3 decimales
  dplyr::mutate_at(vars("p"), funs(round(.,3))) %>% 
  dplyr::mutate(p=ifelse(p<0.001,"<0.001",as.character(p))) %>% 
  dplyr::mutate(CI_95=paste0(sprintf("%3.2f", CI_low),", ", sprintf("%3.2f", CI_high))) %>% 
  dplyr::select(-CI_low, -CI_high, -df_error) 

modelo_final<-
parameters::model_parameters(did.reg) %>% 
  data.frame() %>% 
  #Cambiamos los nombres de las variables
  dplyr::mutate(Parameter= dplyr::case_when(Parameter=="txtime1"~"14 de Febrero hasta el día de obtención",
                                       Parameter=="Primera_100mil"~"Primera dosis por cada 100 mil hab.",
                                       Parameter=="PCR_100mil"~"PCR's por cada 100 mil hab.",
                                       Parameter=="Hosp_UCI_100mil"~"Hospitalizaciones por cada 100 mil hab.",
                                       Parameter=="txtime1:met1"~"Coquimbo y Valparaíso x 14 de Febrero"
                                       )) %>% 
  dplyr::mutate_at(vars("Coefficient","SE","CI_low", "CI_high","t", "df_error"), funs(round(.,2))) %>% 
  #Aproximamos a 3, de manera que los interesados puedan ver el valor p con 3 decimales
  dplyr::mutate_at(vars("p"), funs(round(.,3))) %>% 
  dplyr::mutate(p=ifelse(p<0.001,"<0.001",as.character(p))) %>% 
  dplyr::mutate(CI_95=paste0(sprintf("%3.2f", CI_low),", ", sprintf("%3.2f", CI_high))) %>% 
  dplyr::select(-CI_low, -CI_high, -df_error) 

modelos_combinados<-
modelo_simple %>% 
  #mezclamos las bases de datos sin perder observaciones de ambas, y aquellas filas que calzan, quedarán unidas.
  dplyr::full_join(modelo_final,by="Parameter","_mod_final") %>% 
  dplyr::select(Parameter, Coefficient.x, CI_95.x, p.x, Coefficient.y, CI_95.y, p.y)

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
options(knitr.kable.NA = '')

modelos_combinados %>% 
  #generar tabla en formato html, definiendo título y nombre de columnas
  knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
               caption = paste0("Tabla 1. Coeficientes Diferencia en Diferencias"),
               col.names = c("Variables Independientes","Coeficiente", "IC 95%","Sig","Coeficiente", "IC 95%","Sig"),
align =c('l',rep('c', 101)))%>%
  #Añadimos otra fila como encabezado
  kableExtra::add_header_above(c("","Modelo simple"=3, "Modelo con covariables"=3)) %>% 
  #destaca la quinta fila correspondiente a la intervención para el tratado
  kableExtra::row_spec(2,bold=T,hline_after = T) %>% 
  kableExtra::add_footnote(c("Nota. IC 95%= Intervalos de Confianza al 95%; Datos obtenidos al ", format(Sys.time(), '%d %B, %Y')), notation = "none",threeparttable =T)%>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 13)

• Estamos viendo si las diferencias en casos confirmados entre el periodo post-intervención con el periodo anterior a la intervención (14 feb), es distinto en estas mismas regiones tratadas, respecto las diferencias pre-post de las regiones anteriores.

se <- function(x) sqrt(var(x)/length(x))

Produc_means <- covid19_chile_vacuna_vs_confirmados_std_pcr_uci %>% 
  dplyr::group_by(Region) %>% 
  transmute(y_median = median(conf_100mil),
            y = conf_100mil, 
            fecha = fecha) %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(y_pred = predict(did.reg) + y_median) %>% 
  dplyr::select(-y_median)

#OTRA FORMA DE ESTIMAR LOS PREDICTORES
#https://stackoverflow.com/questions/28547614/prediction-with-plm-method
pmodel_fixed <-pmodel.response(did.reg,model='fixed') 

Produc_means %>% 
  gather(type, value, y, y_pred) %>% 
  dplyr::mutate(tx=factor(ifelse(Region %in% c("Coquimbo", "Valparaíso"),1,0))) %>% 
  dplyr::group_by(tx,type,fecha) %>% 
  dplyr::mutate(median=median(value), se=se(value)) %>%
  dplyr::ungroup() %>% 
  ggplot(aes(x = fecha, y = median, linetype = type))+
  geom_line() +
  geom_ribbon(aes(y= median, x= fecha, ymin = median-se, ymax = median+se), fill = "grey70", color="grey70", alpha=.4)+
  scale_linetype_manual(name="Valores",values=c("solid","dashed"), labels=c("Actual", "Estimado"))+
  facet_wrap(~tx, labeller=as_labeller(c(`0` = "Otras Regiones", `1` = "Coquimbo y\nValparaíso")))+
  geom_vline(xintercept = as.Date("2021-02-14"),linetype="longdash")+
  theme_blank()+
  ggtitle("Prediciones y Medianas por fecha de nuevos casos confirmados")+
  labs(caption="Nota. Área gris representa el error estandar por región;\nLínea vertical= 14 de Febrero",
       x="Fecha",
       y="Mediana Casos Confirmados")

library(gplots)
#defino que se presente un gráfico en una fila, y en 2 columnas #|_|_|
par(mfrow=c(1,2))
# formato función, obtiene una figura con medias de confirmados por 100 mil hab por Región
gplots::plotmeans(conf_100mil ~ Region, main = "Heterogeneidad entre Regiones", ylab="Casos Confirmados (x 100,000)", data = covid19_chile_vacuna_vs_confirmados_std_pcr_uci, n.label=F)
# obtiene una figura con medias de confirmados por 100 mil hab por fechas
gplots::plotmeans(conf_100mil ~ fecha, main = "Heterogeneidad entre Fechas", ylab="", data = covid19_chile_vacuna_vs_confirmados_std_pcr_uci, n.label=F) # no incluyo la etiqueta del eje Y porque ya lo definí en la primera figura.

• A partir del gráfico no es tan claro distinguir si las diferencias se deben al 14 de febrero
• Es claro que las variables de control jugaron un rol confusor aunque falten otros (ej., pasos)
• ¿Cuál es la distribución de los datos?, no lo revisamos
• Las covariables podrían verse afectadas por la intervención (14 feb), en particular en las regiones tratadas
• Adicionalmente, la capacidad de predicción del modelo se estima baja
• La diferencia en cada punto de tiempo no es constante en el periodo pre-intervención
• Sólo dos tratados
• ¿Y la heterogeneidad intraregional? (Comunas distintas)
• Spillover, elementos geográficos no capturados
• Hay patrones temporales que no son capturados (ciclos, tendencias, estaciones, retrasos) (pbgtest())
• Las primeras vacunas siguen un patrón muy distinto
• Volvamos al gráfico anterior (efecto “Valparaíso”, aunque no estuvo contabilizado por 100,000 hab)

#Creamos una función para capturar los intervalos de confianza 
cl_quantile <- function(x, q = c(0.25, 0.75), na.rm = TRUE){
  dat <- data.frame(ymin = quantile(x, probs = q[1], na.rm = na.rm),
                    ymax = quantile(x, probs = q[2], na.rm = na.rm))
  return(dat)
}
  # Defino los parámetros del gráfico en general: no para una capa, sino para todas
ggplot(covid19_chile_vacuna_vs_confirmados_std_pcr_uci, aes(fecha, conf_100mil, color = met)) + 
# Defino las líneas en gris se definieron a partir de datos que excluyen (filtran) aquellas regiones que no son Coquimbo y Valparaíso (met==0), define en la estétitca la fecha en eje x, casos confirmados en el eje y, y define región como grupo (por que si no habrían más datos para cada punto cartográfico igual a N Regiones - Coquimbo y Valpo)
    geom_line(data=covid19_chile_vacuna_vs_confirmados_std_pcr_uci %>% dplyr::filter(met==0), 
              aes(fecha, conf_100mil,group=Region), color="grey80")+
# Hago un resumen de los datos, definiendo que estarán en formato línea y resumirán los datos en la mediana. Fíjese que no tuve que definir el eje x e y ni los datos ni el color de agrupación, ya que los predefiní en la primera línea
    stat_summary(geom = 'line', fun.y = median) + #como línea, la mediana
# Hago otro resumen de los datos, definiendo que estarán en formato de cinta/área y resumirán los datos a partir de la función predefinida al principios. Fíjese que se agregó otro criterio estético: fill (relleno del área), con un valor fijo alpha (transparencia), que tiene un 30% de transparencia ## todo lo que no sea fijo, debe ir en "aes"
    stat_summary(geom = "ribbon", fun.data = cl_quantile, alpha = 0.3, aes(fill=met))+ 
# Agrego una línea vertical en el eje x al 14 de febrero de 2021 #Tiempo de intervención
    geom_vline(xintercept = as.Date("2021-02-14")) + 
  # Defino las etiquetas del color y relleno para la leyenda. #Si asignan los mismos valores, se mezclarán ambas especificaciones: color y relleno del área #Coquimbo y Valparaíso versus el resto
    scale_color_manual(name="Grupo",labels = c("Otras\nRegiones", "Coquimbo y\n Valparaíso"), values=c("red","blue"))+ 
      scale_fill_manual(name="Grupo",labels = c("Otras\nRegiones", "Coquimbo y\n Valparaíso"), values=c("red","blue"))+ 
  # Agrego un tema etilístico para el gráfico
    sjPlot::theme_sjplot2()+
  # Generar etiquetas para el eje x, el eje y y la nota.
  labs(caption="Nota. Área coloreada: rango intecuartil;\nLínea vertical: 14 de febrero 2021;\nLínea intermedia: mediana;\nLíneas en gris: Regiones excluyendo a Coquimbo y Valparaíso",
       y="Casos confirmados (x 100.000 hab.)",
       x="Fecha")

• No son del todo comparables las series (ausencia soporte común)

library(tableone)

#una función que nos permite exportar la tabla en formato kable (html, interactiva)
kableone <- function(x, ...) {
  capture.output(x <- print(x,...))
  knitr::kable(x,
               format= "html", 
               caption = 'Tabla 2. Descriptivos Variables de Interés, antes de 14 de Febrero',
               col.names=c("Resto de\nlas Regiones", "Coquimbo y\nValparaíso","Sig", "Prueba","SMD"), 
               format.args= list(decimal.mark= ".", big.mark= ","))
}

#defino las etiquetas de la base de datos
attr(covid19_chile_vacuna_vs_confirmados_std_pcr_uci$conf_100mil, "label") <- 'Casos confirmados por cada 100 mil hab.'
attr(covid19_chile_vacuna_vs_confirmados_std_pcr_uci$Primera_100mil, "label") <- "Primera dosis por cada 100 mil hab."
attr(covid19_chile_vacuna_vs_confirmados_std_pcr_uci$PCR_100mil, "label") <- "PCR's por cada 100 mil hab."
attr(covid19_chile_vacuna_vs_confirmados_std_pcr_uci$Hosp_UCI_100mil, "label") <- "Hospitalizaciones por cada 100 mil hab."


#tabla comparativa por esrtrato entre tratados (met), se descartan los perdidos (includeNA) (no existen en este caso),
#incorporamos SMD (diferencia de medias estandarizadas)
tab<-
   CreateTableOne(vars = c("conf_100mil","Primera_100mil","PCR_100mil","Hosp_UCI_100mil"), 
                     data= covid19_chile_vacuna_vs_confirmados_std_pcr_uci %>% dplyr::filter(txtime==0), 
                     strata = "met", 
                     includeNA =F,
                      smd=T)

#utilizamos la función predefinida, definiendo que las variables no siguen una distribución normal
kableone(tab, nonnormal= c("conf_100mil","Primera_100mil","PCR_100mil","Hosp_UCI_100mil"),
                       test=T, varLabels=T,noSpaces=T, printToggle=T, dropEqual=F,smd=T) %>% 
  #agrego un formato de las tablas de manera que quede presentable
   kableExtra::kable_styling(bootstrap_options = c("striped", "hover","condensed"),font_size= 12) %>%
  #la primera fila la destacaré con letra negrita porque contiene información relevante
    kableExtra::row_spec(1, bold = T) %>%
  #Agrego una nota
    kableExtra::add_footnote(c("Nota. SMD= Diferencia de Medias Estandarizada"), notation = "none")%>%
  #Si es una tabla muy larga, puedo hacer que sólo abarque un área y lo dejo en formato de tabla deslizable
    kableExtra::scroll_box(width = "100%", height = "200px")